---
title: "Air quality (Covid-19 Research)"
author: "Francesco Bailo"
date: "17/02/2021"
output: html_document
---

# NOTE: Compile from command line 
R -e "Sys.setenv(RSTUDIO_PANDOC='/usr/lib/rstudio-server/bin/pandoc'); rmarkdown::render('README.Rmd',output_file='README.html')"

```{r libraries, message=F, warning=F}
library(dplyr)
library(knitr)
library(sf)
library(ggplot2)
library(viridis)
library(raster)
library(parallel)
library(gridExtra)

knitr::opts_chunk$set(echo = TRUE, message=F, warning=F)
```

```{r}
dates <- seq(as.Date("2020-01-01"), as.Date("2021-01-31"), by = 1)
dates
```

```{r eval = F}

load("../grid/hex_ita_10km.sf.RData")
st_crs(hex_ita_10km.sf) <- 32632
hex_ita_10km.sf <- st_transform(hex_ita_10km.sf, 4326)

fun1 <- function(i, var, hex_ita_10km.sf, dates) {
  require(raster)
  require(sf)
  r <- raster(sprintf('%s.tiff', var), band = i)
  crs(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
  res <- raster::extract(r, hex_ita_10km.sf, fun = mean, weights  = T)
  return(data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                    date = dates[i],
                    var = var,
                    value = res[,1]))
}

cl <- makeCluster(10)

airquality.df <- data.frame(stringsAsFactors = F)

for (var in vars) {
  res.list <- parLapply(cl, 1:length(dates), fun1, var, hex_ita_10km.sf, dates)
  res.df <- bind_rows(res.list)
  airquality.df <- rbind(airquality.df, res.df)
}

stopCluster(cl)

save(airquality.df, file = "airquality.df.RData")

```

# Air quality

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125413.0054502-22845-1-cad9c77f-15c0-47a8-8ce9-10e4d00edcbf-co.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125413.0054502-22845-1-cad9c77f-15c0-47a8-8ce9-10e4d00edcbf-co.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125804.0726807-7813-3-70a95b2c-b7ba-46c1-b5f0-bbc46f157615-pm2p5.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125804.0726807-7813-3-70a95b2c-b7ba-46c1-b5f0-bbc46f157615-pm2p5.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125882.1716464-11197-1-8a9c9eef-8d33-4e22-9e56-a3473de9ed26-so2.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125882.1716464-11197-1-8a9c9eef-8d33-4e22-9e56-a3473de9ed26-so2.tif 

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125644.491702-19159-1-e65235a3-6c25-4200-a8c5-1b02554f0ee1-o3.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125644.491702-19159-1-e65235a3-6c25-4200-a8c5-1b02554f0ee1-o3.tif 

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125496.9399698-5567-1-6822baaf-51d3-4c17-8349-003737620856-no2.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125496.9399698-5567-1-6822baaf-51d3-4c17-8349-003737620856-no2.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125570.481306-21715-6-3b10908e-d9e3-4cb8-b5ee-7b297688ea35-no.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125570.481306-21715-6-3b10908e-d9e3-4cb8-b5ee-7b297688ea35-no.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1614125727.6131365-20235-2-b7ad3e0d-a2c1-4ab6-abc1-591d7fc9bd59-pm10.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1614125727.6131365-20235-2-b7ad3e0d-a2c1-4ab6-abc1-591d7fc9bd59-pm10.tif

```{r eval = F}

load("../grid/hex_ita_10km.sf.RData")

st_crs(hex_ita_10km.sf) <- 32632
require(raster)
require(sf)
require(sp)

dates <- 
  seq(from = as.Date("2020-01-01"), 
      to = as.Date("2021-01-31"),
      by = "day")

pm10.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125727.6131365-20235-2-b7ad3e0d-a2c1-4ab6-abc1-591d7fc9bd59-pm10.tif')
names(pm10.r) <- 
  dates

no.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125570.481306-21715-6-3b10908e-d9e3-4cb8-b5ee-7b297688ea35-no.tif')
names(no.r) <- 
  dates

no2.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125496.9399698-5567-1-6822baaf-51d3-4c17-8349-003737620856-no2.tif')
names(no2.r) <- 
  dates

so2.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125882.1716464-11197-1-8a9c9eef-8d33-4e22-9e56-a3473de9ed26-so2.tif')
names(so2.r) <- 
  dates

pm2p5.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125804.0726807-7813-3-70a95b2c-b7ba-46c1-b5f0-bbc46f157615-pm2p5.tif')
names(pm2p5.r) <- 
  dates

co.r <- 
  brick('adaptor.cams_regional_fc.retrieve-1614125413.0054502-22845-1-cad9c77f-15c0-47a8-8ce9-10e4d00edcbf-co.tif')
names(co.r) <- 
  dates

o3.r <-
  brick('adaptor.cams_regional_fc.retrieve-1614125644.491702-19159-1-e65235a3-6c25-4200-a8c5-1b02554f0ee1-o3.tif')
names(o3.r) <-
  dates

# crs(r) <- 
#   "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
# 
# ita.sp <- 
#   readRDS("gadm36_ITA_0_sp.rds")

hex_ita_10km.sp <- 
  as_Spatial(hex_ita_10km.sf %>% st_transform(4326))

pm10.extr <- 
  raster::extract(pm10.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

pm2p5.extr <- 
  raster::extract(pm2p5.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

co.extr <- 
  raster::extract(co.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

no2.extr <- 
  raster::extract(no2.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

no.extr <- 
  raster::extract(no.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

o3.extr <- 
  raster::extract(o3.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

so2.extr <- 
  raster::extract(so2.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

airquality.df <- data.frame()

for (i in 1:length(dates)) {
  print(i)
  airquality.df <-
    rbind(airquality.df,
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'pm10_conc',
                     value = pm10.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'pm2p5_conc',
                     value = pm2p5.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'co_conc',
                     value = co.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'no2_conc',
                     value = no2.extr[,i]),
          
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'no_conc',
                     value = no.extr[,i]),
          
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'o3_conc',
                     value = o3.extr[,i]),
          
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'so2_conc',
                     value = so2.extr[,i])
    )
}

save(airquality.df, file = "airquality.df.RData")

```

```{r eval = F}
fun2 <- function(i) {
  smpl_var <- sample(unique(airquality.df$var), 1)
  smpl_date <- sample(unique(airquality.df$date), 1)
  smpl_dat <- airquality.df %>% filter(var == smpl_var & date == smpl_date)
  hex_ita_10km.sf$smpl_value <- smpl_dat$value
  return(ggplot(hex_ita_10km.sf) +
           geom_sf(aes(fill = smpl_value), colour = NA) +
           scale_fill_distiller(palette = "Spectral", direction = 1) +
           theme(axis.text = element_blank(),
                 axis.ticks = element_blank()) +
           labs(caption = paste0(smpl_var, " ", smpl_date), fill = NULL) +
           guides(fill = FALSE))
}

plot.list <- lapply(1:25, fun2)
do.call("grid.arrange", c(plot.list, ncol=5))
```


# Temperature & Humidity


gdal_translate -a_srs EPSG:4326 .grib -of Gtiff .tif


Files:
adaptor.mars.internal-1613542488.2219167-23472-25-4599c3c5-f63c-466e-b43b-b1c461bbc8cb-temp.grib
adaptor.mars.internal-1613542518.5545087-6392-3-2ae82254-839a-4f3e-bf3a-f10528d883ac-dewp.grib
adaptor.mars.internal-1613542460.1225214-13302-27-b8efaa2d-2508-45c2-9710-6cbcd0ab550a-prec.grib
adaptor.mars.internal-1613542505.2807376-11678-25-42f58f74-7c7d-4e0c-a249-18c00c2e973d-pres.grib

Commands:
gdal_translate -a_srs EPSG:4326 adaptor.mars.internal-1613542488.2219167-23472-25-4599c3c5-f63c-466e-b43b-b1c461bbc8cb-temp.grib -of Gtiff adaptor.mars.internal-1613542488.2219167-23472-25-4599c3c5-f63c-466e-b43b-b1c461bbc8cb-temp.tif

gdal_translate -a_srs EPSG:4326 adaptor.mars.internal-1613542518.5545087-6392-3-2ae82254-839a-4f3e-bf3a-f10528d883ac-dewp.grib -of Gtiff adaptor.mars.internal-1613542518.5545087-6392-3-2ae82254-839a-4f3e-bf3a-f10528d883ac-dewp.tif

gdal_translate -a_srs EPSG:4326 adaptor.mars.internal-1613542460.1225214-13302-27-b8efaa2d-2508-45c2-9710-6cbcd0ab550a-prec.grib -of Gtiff adaptor.mars.internal-1613542460.1225214-13302-27-b8efaa2d-2508-45c2-9710-6cbcd0ab550a-prec.tif

gdal_translate -a_srs EPSG:4326 adaptor.mars.internal-1613542505.2807376-11678-25-42f58f74-7c7d-4e0c-a249-18c00c2e973d-pres.grib -of Gtiff adaptor.mars.internal-1613542505.2807376-11678-25-42f58f74-7c7d-4e0c-a249-18c00c2e973d-pres.tif

```{r}
load("../grid/hex_ita_10km.sf.RData")
st_crs(hex_ita_10km.sf) <- 32632
require(raster)
require(sf)
require(sp)

dates <- 
  seq(from = as.Date("2020-01-01"), 
      length.out = 409,
      by = "day")

prec.r <- 
  brick('adaptor.mars.internal-1613542460.1225214-13302-27-b8efaa2d-2508-45c2-9710-6cbcd0ab550a-prec.tif')
names(prec.r) <- 
  dates

temp.r <- 
  brick('adaptor.mars.internal-1613542488.2219167-23472-25-4599c3c5-f63c-466e-b43b-b1c461bbc8cb-temp.tif')
names(temp.r) <- 
  dates

pres.r <- 
  brick('adaptor.mars.internal-1613542505.2807376-11678-25-42f58f74-7c7d-4e0c-a249-18c00c2e973d-pres.tif')
names(pres.r) <- 
  dates

dewp.r <- 
  brick('adaptor.mars.internal-1613542518.5545087-6392-3-2ae82254-839a-4f3e-bf3a-f10528d883ac-dewp.tif')
names(dewp.r) <- 
  dates

hex_ita_10km.sp <- 
  as_Spatial(hex_ita_10km.sf %>% st_transform(4326))

temp.extr <- 
  raster::extract(temp.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

prec.extr <- 
  raster::extract(prec.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

pres.extr <- 
  raster::extract(pres.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)

dewp.extr <- 
  raster::extract(dewp.r, hex_ita_10km.sp, 
                  fun = mean, weights  = T)


era5_climate_hex.df <- data.frame()

for (i in 1:length(dates)) {
  print(i)
  era5_climate_hex.df <-
    rbind(era5_climate_hex.df,
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'temp',
                     value = temp.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'prec',
                     value = prec.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'pres',
                     value = pres.extr[,i]),
          data.frame(hex_id = 1:nrow(hex_ita_10km.sf),
                     date = dates[i],
                     var = 'dewp',
                     value = dewp.extr[,i]))
}

save(era5_climate_hex.df, file = "era5_climate_hex.df.RData")


test.sf <- 
  hex_ita_10km.sf

test.sf$value <- dewp.extr[,10]

ggplot(test.sf) + geom_sf(aes(fill = value), colour = NA)

```

# Europe air quality

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1596182248.5360487-22339-5-609faa4c-09e0-4d9a-9dbe-bbe268b3d9a3.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1596182248.5360487-22339-5-609faa4c-09e0-4d9a-9dbe-bbe268b3d9a3.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1596182206.7649863-17613-19-e6d50fd9-23e2-4df2-9246-edd8a5a26a40.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1596182206.7649863-17613-19-e6d50fd9-23e2-4df2-9246-edd8a5a26a40.tif

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1596182217.0965886-2336-10-5d04ee02-99d2-4ec0-a363-c05df4880359.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1596182217.0965886-2336-10-5d04ee02-99d2-4ec0-a363-c05df4880359.tif 

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1596182239.3134377-26085-2-511e1ae7-853a-4e7e-af36-8d90070cd541.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1596182239.3134377-26085-2-511e1ae7-853a-4e7e-af36-8d90070cd541.tif 

gdal_translate -a_srs EPSG:4326 adaptor.cams_regional_fc.retrieve-1596182225.2788684-23525-8-2ed5cc6b-7ba3-4d0c-a153-ff84a9588882.grib -of Gtiff adaptor.cams_regional_fc.retrieve-1596182225.2788684-23525-8-2ed5cc6b-7ba3-4d0c-a153-ff84a9588882.tif 



```{r eval = F}
require(raster)
require(sf)
require(sp)

dates <- seq(from = as.Date("2020-02-11"), to = as.Date("2020-02-23"), by = 1)

europe.sf <- 
  read_sf("../shapefile/Europe/CNTR_BN_60M_2020_4326.shp")

region.sf <- 
  read_sf("../shapefile/Italy/Limiti01012020/Reg01012020/Reg01012020_WGS84.shp")

europe_box.sf <- 
  sf::st_sfc(list(st_polygon(list(matrix(c(-15.1, 29.3, 29.3, -15.1, -15.1,
                           62.1,  62.1,  35.9,  35.9,  62.1), , 2))))) %>%
  st_as_sf( crs = 4326)

no.r <- 
  brick("adaptor.cams_regional_fc.retrieve-1596182206.7649863-17613-19-e6d50fd9-23e2-4df2-9246-edd8a5a26a40.tif")
no.r <-
  crop(no.r, europe_box.sf)

no2.r <- 
  brick("adaptor.cams_regional_fc.retrieve-1596182239.3134377-26085-2-511e1ae7-853a-4e7e-af36-8d90070cd541.tif")
no2.r <-
  crop(no2.r, europe_box.sf)

co.r <-
  brick("adaptor.cams_regional_fc.retrieve-1596182248.5360487-22339-5-609faa4c-09e0-4d9a-9dbe-bbe268b3d9a3.tif")
co.r <-
  crop(co.r, europe_box.sf)

pm2p5.r <-
  brick("adaptor.cams_regional_fc.retrieve-1596182217.0965886-2336-10-5d04ee02-99d2-4ec0-a363-c05df4880359.tif ")
pm2p5.r <-
  crop(pm2p5.r, europe_box.sf)

pm10.r <-
  brick("adaptor.cams_regional_fc.retrieve-1596182225.2788684-23525-8-2ed5cc6b-7ba3-4d0c-a153-ff84a9588882.tif ")
pm10.r <-
  crop(pm10.r, europe_box.sf)

savePlot <- function(r, label) {
  
  these_spdf <- list()
  for (i in 1:length(dates)) {
    these_spdf[[i]] <- as(r[[i]], "SpatialPixelsDataFrame")
    these_spdf[[i]] <- as.data.frame(these_spdf[[i]])
    colnames(these_spdf[[i]]) <- c("value", "x", "y")
  }
  
  plotEur <- function(i, label) {
    ggplot() +
      geom_tile(data = these_spdf[[i]], aes(x=x, y=y, fill=value), alpha=0.8) +
      geom_sf(data = europe.sf, fill = NA, size = .2, alpha = .5, 
              colour = 'white') +
      coord_sf(xlim = c(-15.1, 29.3), ylim = c(35.9, 62.1)) +
      scale_fill_viridis() +
      labs(x = NULL, y = NULL, title = paste0(label," ", format(dates[[i]],"%d %b %Y")), fill = NULL) +
      theme_bw()
  }
  
  these_eur_plots <- lapply(1:12, plotEur, label)
  
  plotNorthIta <- function(i, label) {
    ggplot() +
      geom_tile(data = these_spdf[[i]], aes(x=x, y=y, fill=value), alpha=0.8) +
      geom_sf(data = region.sf %>% st_transform(4326), fill = NA, 
              size = .2, alpha = .5, colour = 'white') +
      coord_sf(ylim = c(43.6, 47.09416), xlim = c(6.612583,14)) +
      scale_fill_viridis() +
      labs(x = NULL, y = NULL, title = paste0(label," ", format(dates[[i]],"%d %b %Y")), fill = NULL) +
      theme_bw()
  }
  
  these_northita_plots <- lapply(1:12, plotNorthIta, label)
  
  ggsave(file = sprintf("../fig/airquality_eur_%s_eur.png", label), width = 18, height = 8,
         arrangeGrob(grobs = these_eur_plots, ncol = 4))
  
  ggsave(file = sprintf("../fig/airquality_northita_%s.png", label), width = 18, height = 8,
         arrangeGrob(grobs = these_northita_plots, ncol = 4))
}

savePlot(pm10.r, label ="PM10")
savePlot(pm2p5.r, label ="PM2.5")
savePlot(co.r, label ="Carbon monoxide")
savePlot(no2.r, label ="Nitrogen dioxide")
savePlot(no.r, label ="Nitrogen monoxide")

```



